#include "cell.h"
#include "config.h"


void Cell::update_T_to_table(   map<double, double> &T_gas,
                                map<double, double> &T_liq)
{
    // status_begin_next.T = status_begin.T;
    // status_mid_next.T = status_mid.T;

    // 此时的 status_begin 已经可以使用
    // 使用此时刻状态，结合比热容 推算出  status_begin_next.T

    // 计算质量
    double m_gas = status_begin.entry.rho_gas * status_begin.alpha_gas * config->cell_volume;
    double m_liq = status_begin.entry.rho_liq * (1 - status_begin.alpha_gas) * config->cell_volume;
    
    // 计算这段时间流失的能量
    double energy_loss = config->heat_coef * config->pipe_perimeter * config->cell_length * 
                        (status_begin.T - config->temp_env) * config->time_step;

    // 计算温度落差
    double delta_T = energy_loss / (m_gas * status_begin.entry.heat_gas + m_liq * status_begin.entry.heat_liq);

    // 更新温度 
    status_begin_next.T = status_begin.T - delta_T;

    // 添加到线性表
    double gas_pos = pos_begin + status_begin.v_gas * config->time_step;
    double liq_pos = pos_begin + status_begin.v_liq * config->time_step;
    T_gas[gas_pos] = status_begin_next.T;
    T_liq[liq_pos] = status_begin_next.T;

    // 打印
    //printf("pos = %lf, T = %lf, T_next = %lf\n", pos_mid, status_mid.T, status_mid_next.T);
}

static double get_map(double pos, map<double, double> &_map){
    auto it = _map.lower_bound(pos);
    if(it == _map.end())
        return (--it)->second;
    if(it == _map.begin())
        return it->second;
    double pos2 = it->first;
    double pos1 = (--it)->first;
    double T1 = it->second;
    double T2 = (++it)->second;
    return T1 + (T2 - T1) / (pos2 - pos1) * (pos - pos1);
}

void Cell::update_T_from_table( map<double, double> &T_gas,
                                map<double, double> &T_liq)
{
    // 计算质量
    double m_gas = status_mid.entry.rho_gas * status_mid.alpha_gas * config->cell_volume;
    double m_liq = status_mid.entry.rho_liq * (1 - status_mid.alpha_gas) * config->cell_volume;
    status_mid_next.T = (status_mid.entry.heat_gas * m_gas * get_map(pos_mid, T_gas) +
                        status_mid.entry.heat_liq * m_liq * get_map(pos_mid, T_liq)) /
                        (status_mid.entry.heat_gas * m_gas + status_mid.entry.heat_liq * m_liq);


    // 同理计算 status_begin_next.T
    m_gas = status_begin.entry.rho_gas * status_begin.alpha_gas * config->cell_volume;
    m_liq = status_begin.entry.rho_liq * (1 - status_begin.alpha_gas) * config->cell_volume;
    status_begin_next.T =   (status_begin.entry.heat_gas * m_gas * get_map(pos_begin, T_gas) +
                            status_begin.entry.heat_liq * m_liq * get_map(pos_begin, T_liq)) /
                            (status_begin.entry.heat_gas * m_gas + status_begin.entry.heat_liq * m_liq);
                                                                                                              
}

void Cell::update_all_to_table( map<double, double> &map_v_gas,
                                map<double, double> &map_v_liq,
                                map<double, double> &map_alpha_gas)
{
    // 添加到线性表
    double gas_pos = pos_begin + status_begin.v_gas * config->time_step;
    double liq_pos = pos_begin + status_begin.v_liq * config->time_step;
    map_v_gas[gas_pos] = status_begin.v_gas;
    map_v_liq[liq_pos] = status_begin.v_liq;
    map_alpha_gas[gas_pos] = status_begin.alpha_gas;
}

void Cell::update_all_from_table(   map<double, double> &map_v_gas,
                                    map<double, double> &map_v_liq,
                                    map<double, double> &map_alpha_gas)
{
    status_begin_next.v_gas = get_map(pos_begin, map_v_gas);
    status_begin_next.v_liq = get_map(pos_begin, map_v_liq);
    status_begin_next.alpha_gas = get_map(pos_begin, map_alpha_gas);

    status_mid_next.v_gas = get_map(pos_mid, map_v_gas);
    status_mid_next.v_liq = get_map(pos_mid, map_v_liq);
    status_mid_next.alpha_gas = get_map(pos_mid, map_alpha_gas);
}

void Cell::update_entry(Status& status)
{
    PT_Entry &entry = status.entry;
    if(entry.P == status.P && entry.T == status.T)
        return;
    entry = pt_table->find(status.P, status.T);
}

void Cell::update_entrance_alpha(Status &status)
{
    update_entry(status);

    // 若仅提供了总输量，则使用平衡气体质量分数分配入口的油输量和气输量
    if(config->vm_total_only)
    {
        config->vm_gas = config->vm_total * status.entry.rs;
        config->vm_liq = config->vm_total * (1 - status.entry.rs);
    }

    double vv_gas = config->vm_gas / status.entry.rho_gas;  // 入口气体体积流量 m^3/s
    double vv_liq = config->vm_liq / status.entry.rho_liq;  // 入口液体体积流量 m^3/s
    double alpha_gas = vv_gas / (vv_gas + vv_liq);  // 初始气体体积分数
    double alpha_liq = vv_liq / (vv_gas + vv_liq);  // 初始液体体积分数
    
    double v_gas = 0;  // 入口气体速度 m/s
    if(alpha_gas > 0)
        v_gas = vv_gas / (config->pipe_area * alpha_gas);
    double v_liq = 0;  // 入口液体速度 m/s
    if(alpha_liq > 0)
        v_liq = vv_liq / (config->pipe_area * alpha_liq);  

    status.v_gas = v_gas;
    status.v_liq = v_liq;
    status.alpha_gas = alpha_gas;
}

double Cell::get_gas_mass_residual()
{
    //等式左侧 量纲 kg/(m^3·s)
    double left_side = (status_mid_next.alpha_gas * status_mid_next.entry.rho_gas -
                        status_mid.alpha_gas * status_mid.entry.rho_gas) /
                        config->time_step;

    //等式右侧 量纲 kg/(m^3·s)
    double right_side = - (right->status_begin_next.alpha_gas * right->status_begin_next.entry.rho_gas * right->status_begin_next.v_gas -
                           status_begin_next.alpha_gas * status_begin_next.entry.rho_gas * status_begin_next.v_gas) /
                           config->cell_length;
          
    //相间传质 量纲 kg/(m^3·s)
    psi_gas = 0;
    psi_gas += (right->status_begin_next.entry.rho_gas * right->status_begin_next.v_gas * right->status_begin_next.alpha_gas +
                right->status_begin_next.entry.rho_liq * right->status_begin_next.v_liq * (1-right->status_begin_next.alpha_gas))
                * right->status_begin_next.entry.rs;
    psi_gas -= (status_begin_next.entry.rho_gas * status_begin_next.v_gas * status_begin_next.alpha_gas +
                status_begin_next.entry.rho_liq * status_begin_next.v_liq * (1-status_begin_next.alpha_gas))
                * status_begin_next.entry.rs;
    psi_gas /= config->cell_length;

    // // 当前气体质量分数
    // double m_gas = status_mid_next.alpha_gas * status_mid_next.entry.rho_gas * config->cell_volume;
    // double m_liq = (1-status_mid_next.alpha_gas) * status_mid_next.entry.rho_liq * config->cell_volume;
    // double m_alpha_gas = m_gas / (m_gas + m_liq);
    // // 基于平衡质量分数的动力学模型
    // double tao = 10;  // 时间常数
    // psi_gas = (status_mid_next.entry.rs - m_alpha_gas) / tao * status_mid_next.entry.rho_liq;

    // 计算等式误差
    double residal = left_side - right_side - psi_gas;
    return residal;
}

double Cell::get_liq_mass_residual()
{
    //等式左侧 量纲 kg/(m^3·s)
    double left_side = ((1-status_mid_next.alpha_gas) * status_mid_next.entry.rho_liq -
                        (1-status_mid.alpha_gas) * status_mid.entry.rho_liq) /
                        config->time_step;

    //等式右侧 量纲 kg/(m^3·s)
    double right_side = - ((1-right->status_begin_next.alpha_gas) * right->status_begin_next.entry.rho_liq * right->status_begin_next.v_liq -
                           (1-status_begin_next.alpha_gas) * status_begin_next.entry.rho_liq * status_begin_next.v_liq) /
                           config->cell_length;

    // 已经计算过了 psi_gas

    // 计算等式误差
    double residal = left_side - right_side + psi_gas;
    return residal;
}

double Cell::get_gas_momentum_residual()
{
    //等式左侧 量纲 kg/(m^2·s^2)
    double left_side = 0;
    left_side += status_mid_next.alpha_gas * status_mid_next.entry.rho_gas * status_mid_next.v_gas;
    left_side -= status_mid.alpha_gas * status_mid.entry.rho_gas * status_mid.v_gas;
    left_side /= config->time_step;

    //等式右侧1 量纲 kg/(m^2·s^2)
    double right_side1 = 0;
    right_side1 += right->status_begin_next.P - status_begin_next.P;
    right_side1 /= config->cell_length;
    right_side1 *= -status_mid_next.alpha_gas;

    //等式右侧2 量纲 kg/(m^2·s^2)
    double right_side2 = 0;
    right_side2 += right->status_begin_next.alpha_gas * right->status_begin_next.entry.rho_gas * right->status_begin_next.v_gas * right->status_begin_next.v_gas;
    right_side2 -= status_begin_next.alpha_gas * status_begin_next.entry.rho_gas * status_begin_next.v_gas * status_begin_next.v_gas;
    right_side2 /= -config->cell_length;

    // 近似气液比例即为周长占比
    double circum_gas = config->pipe_perimeter * status_mid_next.alpha_gas;
    double circum_liq = config->pipe_perimeter * (1-status_mid_next.alpha_gas);
    double circum_gas_liq = circum_gas < circum_liq ? circum_gas : circum_liq; // 近似为气液交界处周长

    // 气体受到管壁的剪切力
    double Re = status_mid_next.entry.rho_gas * abs(status_mid_next.v_gas) * (2 * config->pipe_radius) / status_mid_next.entry.visc_gas;  // 雷诺数
    double f;  // 摩擦系数
    // 根据雷诺数判断流动状态并计算摩擦系数
    if (Re == 0)
        f = 0;
    else if (Re < 2300) // 层流
        f = 64 / Re;
    else    // 湍流，使用布拉休斯公式
        f = 0.316 * pow(Re, -0.25);
    // 计算壁面剪切力 (Kg/m*s^2)
    double tau_wall = 0.5 * f * status_mid_next.entry.rho_gas * status_mid_next.v_gas * abs(status_mid_next.v_gas);

    //等式右侧3 量纲 kg/(m^2·s^2)
    double right_side3 = -tau_wall * circum_gas / config->pipe_area;

    // 计算气体受液体的剪切力
    double v_rel = status_mid_next.v_gas - status_mid_next.v_liq;  // 相对速度
    double rho_m = status_mid_next.alpha_gas * status_mid_next.entry.rho_gas + (1-status_mid_next.alpha_gas) * status_mid_next.entry.rho_liq;  // 混合密度
    double visc_m = status_mid_next.alpha_gas * status_mid_next.entry.visc_gas + (1-status_mid_next.alpha_gas) * status_mid_next.entry.visc_liq;  // 混合粘度
    // 计算界面摩擦系数
    double Re_i = rho_m * abs(v_rel) * config->pipe_diameter / visc_m;  // 界面雷诺数
    double f_i;  // 界面摩擦系数
    if (Re_i == 0)
        f_i = 0;
    else if (Re_i < 2300) {
        f_i = 64 / Re_i;  // 层流
    } else {
        f_i = 0.316 * pow(Re_i, -0.25);  // 湍流
    }
    // 计算界面剪切力
    double tau_interface = 0.5 * f_i * rho_m * v_rel * abs(v_rel);

    // 等式右侧4 量纲 kg/(m^2·s^2)
    double right_side4 = -tau_interface * circum_gas_liq / config->pipe_area;

    // 等式右侧5 量纲 kg/(m^2·s^2)
    double va = psi_gas > 0 ? status_mid_next.v_liq : status_mid_next.v_gas;
    double right_side5 = psi_gas * va;

    // 等式右侧6 量纲 kg/(m^2·s^2)
    double right_side6 = -status_mid_next.alpha_gas * status_mid_next.entry.rho_gas * config->g * sin_theta;

    // 计算等式误差
    double residal = left_side - right_side1 - right_side2 - right_side3 - right_side4 - right_side5;
    return residal;
}

double Cell::get_liq_momentum_residual()
{
    //等式左侧 量纲 kg/(m^2·s^2)
    double left_side = 0;
    left_side += (1-status_mid_next.alpha_gas) * status_mid_next.entry.rho_liq * status_mid_next.v_liq;
    left_side -= (1-status_mid.alpha_gas) * status_mid.entry.rho_liq * status_mid.v_liq;
    left_side /= config->time_step;

    //等式右侧1 量纲 kg/(m^2·s^2)
    double right_side1 = 0;
    right_side1 += right->status_begin_next.P - status_begin_next.P;
    right_side1 /= config->cell_length;
    right_side1 *= -(1-status_mid_next.alpha_gas);

    //等式右侧2 量纲 kg/(m^2·s^2)
    double right_side2 = 0;
    right_side2 += (1-right->status_begin_next.alpha_gas) * right->status_begin_next.entry.rho_liq * right->status_begin_next.v_liq * right->status_begin_next.v_liq;
    right_side2 -= (1-status_begin_next.alpha_gas) * status_begin_next.entry.rho_liq * status_begin_next.v_liq * status_begin_next.v_liq;
    right_side2 /= -config->cell_length;

    // 近似气液比例即为周长占比
    double circum_gas = config->pipe_perimeter * status_mid_next.alpha_gas;
    double circum_liq = config->pipe_perimeter * (1-status_mid_next.alpha_gas);
    double circum_gas_liq = circum_gas < circum_liq ? circum_gas : circum_liq; // 近似为气液交界处周长

    // 液体受到管壁的剪切力
    double Re = status_mid_next.entry.rho_liq * abs(status_mid_next.v_liq) * (2 * config->pipe_radius) / status_mid_next.entry.visc_liq;  // 雷诺数
    double f;  // 摩擦系数
    // 根据雷诺数判断流动状态并计算摩擦系数
    if (Re == 0)
        f = 0;
    else if (Re < 2300) // 层流
        f = 64 / Re;
    else    // 湍流，使用布拉休斯公式
        f = 0.316 * pow(Re, -0.25);
    // 计算壁面剪切力 (Kg/m*s^2)
    double tau_wall = 0.5 * f * status_mid_next.entry.rho_liq * status_mid_next.v_liq * abs(status_mid_next.v_liq);

    //等式右侧3 量纲 kg/(m^2·s^2)
    double right_side3 = -tau_wall * circum_liq / config->pipe_area;

    // 计算液体受气体的剪切力
    double v_rel = status_mid_next.v_liq - status_mid_next.v_gas;  // 相对速度
    double rho_m = status_mid_next.alpha_gas * status_mid_next.entry.rho_gas + (1-status_mid_next.alpha_gas) * status_mid_next.entry.rho_liq;  // 混合密度
    double visc_m = status_mid_next.alpha_gas * status_mid_next.entry.visc_gas + (1-status_mid_next.alpha_gas) * status_mid_next.entry.visc_liq;  // 混合粘度
    // 计算界面摩擦系数
    double Re_i = rho_m * abs(v_rel) * config->pipe_diameter / visc_m;  // 界面雷诺数
    double f_i;  // 界面摩擦系数
    if (Re_i == 0)
        f_i = 0;
    else if (Re_i < 2300) {
        f_i = 64 / Re_i;
    } else {
        f_i = 0.316 * pow(Re_i, -0.25);
    }
    // 计算界面剪切力
    double tau_interface = 0.5 * f_i * rho_m * v_rel * abs(v_rel);

    // 等式右侧4 量纲 kg/(m^2·s^2)
    double right_side4 = -tau_interface * circum_gas_liq / config->pipe_area;

    // 等式右侧5 量纲 kg/(m^2·s^2)
    double va = psi_gas > 0 ? status_mid_next.v_liq : status_mid_next.v_gas;
    double right_side5 = -psi_gas * va;

    // 等式右侧6 量纲 kg/(m^2·s^2)
    double right_side6 = -(1-status_mid_next.alpha_gas) * status_mid_next.entry.rho_liq * config->g * sin_theta;

    // 等式右侧7 量纲 kg/(m^2·s^2)
    double right_side7 = -(1-status_mid_next.alpha_gas) * config->pipe_diameter *
                        (status_mid_next.entry.rho_liq - status_mid_next.entry.rho_gas) 
                        * config->g * sin_theta *
                        ((1-right->status_begin_next.alpha_gas)-(1-status_begin_next.alpha_gas)) /
                        config->cell_length;

    // 计算等式误差
    double residal = left_side -    right_side1 - right_side2 - right_side3 - 
                                    right_side4 - right_side5 - right_side6 - 
                                    right_side7;
    return residal;
}


double Cell::get_gas_inlet_residual()
{
    // 入口气体定流的残差即为 应有的入口气体速度与实际入口气体速度的差值
    double vv_gas = config->vm_gas / status_begin_next.entry.rho_gas;  // 应有的入口气体体积流量 m^3/s
    double v_gas = vv_gas / (config->pipe_area * status_begin_next.alpha_gas);  // 应有的入口气体速度 m/s
    if(config->vm_gas == 0)  // 若入口无气体流量，则速度为0
        v_gas = 0;
    double residal = status_begin_next.v_gas - v_gas;
    return residal;
}

double Cell::get_liq_inlet_residual()
{
    // 入口液体定流的残差即为 应有的入口液体速度与实际入口液体速度的差值
    double vv_liq = config->vm_liq / status_begin_next.entry.rho_liq;  // 应有的入口液体体积流量 m^3/s
    double v_liq = vv_liq / (config->pipe_area * (1-status_begin_next.alpha_gas));  // 应有的入口液体速度 m/s
    if(config->vm_liq == 0)  // 若入口无液体流量，则速度为0
        v_liq = 0;
    double residal = status_begin_next.v_liq - v_liq;
    return residal;
}

double Cell::get_outlet_P_residual()
{
    // 出口定压的残差即为 出口压强与右侧网格的压强差值
    double residal = config->P_exit - right->status_begin_next.P;
    return residal;
}

// 将下一时刻状态应用到当前状态
void Cell::apply_next_status()
{
    status_begin = status_begin_next;
    status_mid = status_mid_next;
}

// 将当前状态应用到下一时刻状态
void Cell::update_next_status()
{
    status_begin_next = status_begin;
    status_mid_next = status_mid;
}

void Cell::update_status_mid_next()
{
    // 若右侧为空，则不更新
    if(right == nullptr)
        return;

    // 若右侧不为空，则更新中间状态为左侧状态与右侧状态的平均值
    status_mid_next.P = 0.5 * (status_begin_next.P + right->status_begin_next.P);
    status_mid_next.T = 0.5 * (status_begin_next.T + right->status_begin_next.T);
    status_mid_next.v_gas = 0.5 * (status_begin_next.v_gas + right->status_begin_next.v_gas);
    status_mid_next.v_liq = 0.5 * (status_begin_next.v_liq + right->status_begin_next.v_liq);
    status_mid_next.alpha_gas = 0.5 * (status_begin_next.alpha_gas + right->status_begin_next.alpha_gas);
    update_entry(status_mid_next);
}